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Abstract 

^  The  transient  heat  conduction  equation  with  a  nonlinear 
conductivity  is  solved  using  three  finite  difference  schemes 
and  two  finite  element  schemes.  The  schemes  are  compared 
for  speed  and  accuracy  against  an  analytic  solution.  The 
finite  element  methods  are  shown  to  be  the  most  accurate  and 
one,  a  lumped  mass  method,  is  comparable  in  speed  to  the 
finite  difference  methods.  The  lumped  mass  finite  element 
method  is  shown  to  be  the  clear  choice  from  among  the 
methods  tested  in  solving  the  test  problem. 


A  COMPARISON  OF  THE  FINITE  DIFFERENCE 


AND  FINITE  ELEMENT  METHODS  FOR  THE  SOLUTION 
OF  THE  TRANSIENT  HEAT  CONDUCTION  EQUATION 
WITH  TEMPERATURE- DEPENDENT  CONDUCTIVITY 

I .  Introduction 

Like  any  differential  equation  that  corresponds  to  a 
real  problem,  the  heat  conduction  equation  only  approximates 
the  physical  phenomena  of  heat  conduction.  The  equation 
itself,  however,  is  usually  further  simplified  by  assuming 
that  the  specific  heat  and  conductivity  of  the  material  are 
constant.  This  assumption  yields  linear  equations  which 
are  much  easier  to  solve  than  the  nonlinear  equations  that 
result  from  the  unsimplified  heat  equation.  In  many  cases, 
the  loss  of  accuracy  is  small  and  is  completely  overshadowed 
by  the  advantage  of  increased  efficiency.  Unfortunately, 
however,  this  is  not  always  true. 

An  example  of  a  highly  nonlinear  problem  is  that  of  a 
fuel  pin  in  a  nuclear  power  reactor.  At  the  center  of  the 
fuel  pin,  the  temperature  is  on  the  order  of  4000°  F  under 
normal  operating  conditions,  and  the  exterior  is  at  around 
600°  F  (Ref  11:23-4).  This  large  temperature  drop  results 
in  large  variation  in  conductivity  over  a  distance  of  about 
half  an  inch.  This  problem  cannot  be  accurately  treated  by 
assuming  a  constant  conductivity,  and  thus  the  nonlinear 
heat  equation  must  be  solved  (Ref  4:104). 
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The  common  techniques  used  by  engineers  to  solve  boundary 
value  problems  in  differential  equations  on  a  computer  are 
finite  differences  and  finite  elements.  These  methods  can  be 
applied  to  nonlinear  heat  conduction  and  one  of  the  two  would 
probably  be  the  choice  for  most  engineers  in  solving  the 
equation  on  a  computer.  For  this  reason,  it  is  useful  to 
know  if  one  method  is  significantly  superior  to  the  other  in 
terms  of  accuracy  and  efficiency. 

There  are  various  ways  to  compare  the  performance  of 
different  numerical  techniques  that  are  applicable  to  a 
problem.  The  most  rigorous  and  satisfying  is  a  mathematical 
treatment  that  yields  figures  of  merit  that  can  be  compared 
among  the  methods.  This  is  not  always  possible  without 
including  qualifications  that  limit  the  scope  of  the  results. 
Another  approach,  which  is  the  one  used  by  the  author,  is 
to  apply  the  numerical  methods  to  a  problem  with  a  known 
solution.  By  measuring  the  effectiveness  of  the  methods  in 
solving  the  problem,  one  can  get  an  idea  if  one  method  is 
significantly  better  than  the  rest.  The  obvious  problem 
with  the  latter  approach  is  that  any  type  of  generalization 
to  the  larger  class  of  problems  is  difficult.  For  this 
reason,  the  results  must  be  viewed  cautiously. 

The  objective  of  this  work  is  to  compare  different 
variations  of  the  finite  difference  and  finite  element 
methods  by  having  them  solve  a  nonlinear  heat  transfer 
problem  with  a  known  solution.  The  methods  will  be  compared 
primarily  on  the  basis  of  accuracy  and  speed. 
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Because  of  the  nature  of  the  test  problem,  a  one¬ 
dimensional  form  of  the  nonlinear  heat  equatijn  will  be 
used  that  hac  a  temperature  dependent  conductivity,  but  a 
constant  specific  heat.  The  comparison  is  thus  limited  to 
one-dimensional  finite  differences  and  finite  elements. 
There  are  also  many  variations  of  the  finite  element  method 
that  can  be  used,  but  in  this  work  those  variations  are 
limited  to  CO  basis  functions  with  Galerkin’s  method.  For 
the  finite  difference  method,  three  differencing  schemes 
are  compared. 
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I I .  Theory 

The  theory  behind  this  work  can  be  divided  into  two 
categories.  The  first  is  a  general  background  of  the  finite 
difference  and  finite  element  techniques  and  the  methods 
used  by  the  author.  In  this  case,  an  understanding  of  the 
techniques  is  presumed  of  the  reader  and  the  presentation 
is  an  outline  with  those  points  peculiar  to  this  work  stressed. 

The  second  general  area  is  the  nonlinear  heat  equation 
with  an  anlaytic  solution  that  will  be  used  as  the  test 
problem.  The  derivation  of  the  solution  is  not  easy  and  will 
not  be  included  in  this  work.  Significant  points  of  the  problem 
are  stressed.  Because  the  evaluation  of  the  analytic  solution 
is  not  trivial  and  must  be  done  numerically,  it  is  included 
in  the  Appendix. 

The  Time-Dependent  Finite  Difference  and  Finite  Element  Methods 

The  finite  difference  and  finite  element  methods  are 
both  used  to  solve  the  transient  nonlinear  heat  conduction 
problem.  When  applied  to  this  problem,  they  both  yield  a 
set  of  equations  that  can  be  written  in  matrix  notation  as 

[A3  =  [B(e)3  0  (1) 


where  F  is  a  vector  of  nodal  temperatures  and  t  is 
time  (Ref  6:135).  Note  that  the  matrix  [B(0)3  is  dependent 
on  9  so  the  equation  is  nonlinear.  These  equations  can  be 
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treated  in  the  same  manner  for  both  finite  differences  and 
finite  elements.  For  this  reason,  the  spatial  discretization 
process  for  each  method  will  be  presented  first,  which  yields 
Eq  (1),  and  then  methods  of  handling  Eq  (1)  will  be  given. 

Finite  Differences 

The  time-dependent  nonlinear  heat  conduction 
equation  to  be  solved  is  given  by 

ST  '  4  <25 

where  0  is  temperature,  t  is  time,  k(0)  is  the  temper¬ 
ature  dependent  conductivity,  and  x  is  the  spatial  coordinate 
It  is  assumed  that  the  specific  heat  and  density  are  unity. 

The  use  of  standard  central  difference  techniques  on  the 
linear  heat  equation  is  straightforward  and  is  normally  done 
in  only  one  way.  By  rewriting  the  nonlinear  heat  equation  in 
other  mathematically  equivalent  forms,  however,  it  is  possible 
to  come  up  with  at  least  two  different  variations  that  appear 
to  be  equally  accurate.  Other  forms  of  the  nonlinear  heat 
equation  that  can  be  used  are 
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3k (0)  30 


3  2  0 


The  finite  difference  equations  for  both  of  these  variations 
will  be  presented. 

Central  difference  formulas  for  the  spatial  domain  are 
well  known,  but  will  be  summarized  here.  To  approximate 
the  first  derivative  of  a  function  at  a  node  i,  one  can  use 


ei+i  •  ei 
Ax-; 


Axi-1 


6i  -  9i-l 
Axj-i 


Ax^ 


Axi+1  +  Axi 


(5) 


where  Ax.  =  x .  . .  -  x. 

1  l+l  l 


and 


Ax .  .  =  x .  -  x .  , 
l-l  l  l-l 


The 


approximation  of  the  second  derivative  is  given  by 


9i+i  -  9i 

AXj 


9i  -  9i-l 
Axl-1 


h CAxi_i  +  AxA) 


(6) 


These  equations  are  applicable  when  the  intervals  are  not 
equally  sized,  i.e.,  Ax^  Ax^ 

To  derive  the  finite  difference  equations,  one  must 
substitute  Eqs  (5)  and  (6)  into  the  different  forms  of  the 
heat  equation.  For  some  node  i  and  Eq  (4), 
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+  keep 

where  it  is  assumed  that  the  analytic  form  for  is 

known.  This  equation  is  nonlinear  because  of  the  first  term 
on  the  right  hand  side. 

For  Eq  (3),  one  gets 


These  two  forms  will  be  referred  to  as  Finite  Difference 
Schemes  1  and  2,  respectively.  It  is  not  evident  a  priori 
that  one  form  is  better  than  the  other. 
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Another  set  of  difference  equations  is  presented  by 
Myers  (Ref  9:298).  These  equations  are  derived  from  an 
energy  balance  based  on  the  heat  flux  given  by 


qi-l,i 


qi,i+l 


k.  +  k. 
l-l  1 


k.  +  k.  . 
i  i+I 

- 2 - 


ei  -  6i-i 
Axi-i 


ei+i  -  et 
Ax^ 


(9) 


(10) 


where  q  is  the  heat  flux,  k  is  the  conductivity,  Ax  is 
the  spatial  increment.  The  resulting  difference  equation  is 
given  by 


,  A  k .  i  +  k  • 

df '  _  ,  l-l  1-.  0 

at yi  "  ( - 2 - 3  i-1 


+ 


+ 

~2 


0i  +  l 


*  ki+i. 


e 

i 


(11 


and  will  be  referred  to  as  Finite  Difference  Scheme  3. 

All  of  the  schemes  presented  above  can  be  written  in 
general  as 

[a:  e  =  CB(e)Je  (i 

where  [AD  is  the  identity  matrix.  This  is  an  initial  value 
problem  in  time,  and  solution  methods  which  are  applicable  to 
both  the  finite  differences  and  finite  elements  will  be  pre¬ 


sented  later. 


Finite  Elements 


The  foundation  of  the  finite  element  method  is  to 
find  a  linear  combination  of  basis  functions  that  approximate 
the  solution  to  a  differential  equation.  The  finite  element 
method  differs  from  similar  techniques  because  the  basis 
functions  are  locally  defined.  This  results  in  relatively 
sparse  matrices  which  are  suitable  for  computer  calculations. 

The  basis  functions  are  selected  with  reference  to 
elements,  which  are  subdivisions  of  the  domain  of  the  prob¬ 
lem.  In  a  one-dimensional  domain,  the  elements  are  simply 
intervals,  as  shown  in  Figure  1.  The  basis  functions  used 
in  this  work  have  continuity  Cc  ,  which  means  that  globally 
they  are  zeroeth-order  differentiable.  They  are  also  piece- 
wise  linear.  The  basis  functions  are  shown  in  Figure  1  and 
are  given  in  an  interval  by 


Na  (x) 


Sx.  s  X  s 


Ni+l(x)  = 


X  . 
1 


xi+l  ■  xi 


X  . 


i+1 


/ 

These  basis  functions  are  sufficiently  complete  for  one¬ 
dimensional  heat  transfer  (Ref  14:64).  Note  that  these 
interpolation  functions  are  not  time-dependent.  If  one 
defines  Tl(t)  and  T2(t)  as  the  temperature  at  the  end¬ 
points  of  an  interval,  then  the  temperature  at  any  point 
x  in  the  interval  is  given  by 


(12) 
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Figure  1.  Local  Basis  Functions  in  an  Interval  or  Element 


T(x)  =  [NjCx)  N.+1(x)] 

Thus  far,  the  basis  functions  have  been  defined  within 
an  interval,  and  hence  on  a  local  basis  only.  One  can 
associate  a  global  interpolation  function  with  every  node. 
The  global  function  is  given  over  the  domain  by 


where  the  basis  function  is  zero  everywhere  except  around 
the  node.  The  approximate  solution  over  the  domain  is 
then  given  by 

e(x,t)  :  N.(x)  0i(t)  =  [N.(x)  9i(t)  (15 

i 

where  6^  are  the  nodal  unknowns  and  summation  is  implied 
over  identical  indices. 

To  determine  the  approximate  solution,  the  nodal  unknowns 
0^  must  be  found.  Substituting  the  approximations  into  the 
heat  conduction  equation,  one  gets 

36.  (t)  3N.  (x)  \ 

"i*”)  -it—  -  sr (kfe)  -t —  vty  *  '16 

where  the  term  R(x,t)  is  the  error,  or  residual,  that 
results  from  introducing  the  approximation.  The  method  of 
weighted  residuals  now  requires  the  selection  of  a  weighting 


V 


a 


I 


function  Wj,  .  In  Galerkin's  method,  a  special  case  of 
the  method  of  weighted  residuals,  the  weighting  functions 
are  just  the  basis  functions,  or  .  Multiplying 

both  sides  of  the  differential  equation  by  the  weighting 
function  and  integrating  over  the  domain,  D  , 


/D  N.  (x)  Nj  (x)  dx 


36.  (t) 
3t 


-  /D  Nj  M 


3x 


(k(6) 


3Ni(x) 

3x 


e.(t) 


) 


dx 


+  JD  R(x,t)  N.  (x)  dx 


(17) 


The  residual  R(x,t)  is  made  orthogonal  to  the  weighting 
functions,  i.e., 


/D  R(x,t)  Nj  (x)  dx  = 


(18) 


and  thus 


3  8- 


8N. 


h>  Ni  Nj  dx  TT  *  Id  Nj  3x  ST  d* 


(19) 


Integrating  by  parts, 
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/p.  N.  N.  dx 
'D  1  j 


3N.  3N. 

0ict’  d* 


3N . 

+  Nj  (x)  --L  kce)  e.(t)  (20) 

This  equation  is  nonlinear  because  the  conductivity  is 
temperature  dependent  on  the  right  hand  side. 

The  integrals  in  Eq  (20)  are  all  known  quantities  and 
can  be  represented  by  matrices  as 


CAD  =  /D  Ni  N^.  dx 

3N.  3N .  3N. 

CBfe)]  -  -/„  Me)  dx  *  N^kce)  (21) 


Rewriting  Eq  (20)  in  matrix  form. 


CAD  0  =  CB(6)D  e 


(1) 


Using  terminology  from  structures,  Eq  (20)  represents 
the  distributed  mass  form  of  the  finite  element  method,  or 
perhaps  more  appropriately  here,  "distributed  heat  capacity." 
This  nomenclature  refers  to  the  averaging  process  and  is 
intimately  related  to  the  choice  of  weighting  functions 
(Ref  14:535).  Another  finite  element  implementation,  the 
lumped  mass  form,  can  also  be  used.  One  method  of  deriving 
lumped  mass  equations  is  to  add  the  off-diagonal  elements 
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on  each  row  of  CAD  to  the  diagonal  elements  to  form  a 
diagonal  matrix  (Ref  14:535).  This  method  was  used  in  this 
work.  In  the  context  of  this  work,  these  two  methods  will 
be  viewed  merely  as  two  different  finite  element  implemen¬ 
tations  and  the  physical  implications  will  not  be  examined. 

The  two  methods  will  be  referred  to  as  the  distributed  mass 
and  lumped  mass  finite  element  methods  in  spite  of  the  implied 
reference  to  structures. 

The  equations  for  both  of  the  finite  element  methods  are 
similar  in  form  to  the  results  of  the  spatial  discretization 
using  finite  differences.  Methods  to  solve  the  system  are 
presented  in  the  next  section. 

Timestepping  Techniques 

Both  the  finite  difference  and  finite  element 
methods  yield  a  set  of  equations  that  can  be  represented 
in  matrix  notation  as 

CA]  9  =  [B(0)3  0  (1 

For  finite  differences,  CA]  is  the  identity  matrix.  It 
has  been  shown  how  these  matrices  are  found  for  both  methods. 
The  next  step  is  to  use  an  initial  value  solution  technique 
to  solve  this  set  of  equations. 

The  standard  algorithm  that  is  generally  used  for 


finite  difference  and  finite  element  time- stepping  is  given 
by  (Ref  6:2) 
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Tn  :  e  (tn)  (22g) 

Un  :  0  (tn)  (22h) 

For  0  s  a  s  1  ,  different  schemes  result.  When  a  =  0  , 

one  has  a  purely  explicit  or  Euler's  method.  When  a  =  1  , 

a  purely  implicit  method  results.  Both  of  these  methods 
have  accuracy  of  the  first  order.  A  second  order  scheme, 
called  the  Crank-Nicholson  method,  results  when  a  =  0.5 
This  method  is  generally  the  most  accurate  of  this  family 
of  schemes  and  was  the  one  used  by  the  author. 


t: 
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Other  considerations  are  important  when  one  selects  an 
initial  value  problem  solver  for  use  with  finite  differences 
or  finite  elements.  A  major  factor  is  the  stability  of 
the  method.  A  method  is  not  stable  if  errors  begin  to  grow 
exponentially.  It  can  be  shown  that  the  Cr  mk- Ni cholson 
method  is  always  stable  for  linear  heat  conduction  (Ref  2:516). 
The  only  method  of  the  family  shown  above  that  is  unconditionally 
stable  for  nonlinear  heat  conduction  is  the  purely  implicit 
method.  This  is  unfortunate  because  the  purely  implicit 
method  only  has  first  order  accuracy  (Ref  6:137).  The  methods 
that  are  in  the  family  represented  by  Eq  (22)  and  are  only 
conditionally  stable  can  be  made  stable  by  limitations  on 
the  time-step  length  for  linear  problems.  Richtmeyer  gives 
formulas  the  limitation  of  the  time  step  length  for  the 
linear  case  (Ref  10:187+).  Hughes  presents  an  analysis  of 
time  step  limitations  for  nonlinear  heat  transfer,  but  does 
not  give  explicit  formulae  (Ref  6:4).  The  author  selected  a 
maximum  time  step  length  using  the  formula  given  by  Richtmeyer 
for  the  linear  case  and  a  Crank-Nicholson  scheme  given  by 
(Ref  10:187) 

“  -  *frr  (23) 

where  At  is  the  time  step  length.  Ax  is  the  spatial 
increment,  |kj  is  the  maximum  conductivity  over  the  domain, 
and  y  is  a  parameter  such  that  0  <  y  <  1  .  For 


y  :  o.i 


,  stability  was  achieved  in  most  cases. 


The  Test  Problem 

Luikov  presents  three  analytic  solutions  to  the 
nonlinear  heat  equation 


30 

3t 


(2) 


where  6  is  temperature,  k(0)  is  the  temperature  dependent 
conductivity,  x  is  the  independent  spatial  variable,  and 
t  is  the  independent  time  variable  (Ref  7:499+).  In  this 
form,  it  is  assumed  that  the  specific  heat  and  density  are 
constant  or  unity.  Luikov ’s  solutions  are  for  semi- inf inite 
domains,  specific  functions  of  conductivity  with  temperature, 
and  specific  boundary  conditions. 

The  nonlinear  heat  equation  and  boundary  conditions  that 
are  used  as  a  benchmark  in  this  paper  are  from  Luikov’s 
first  problem,  which  was  originally  presented  by  Fujita 
(Ref  7:499;  5),  and  is  given  by 
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0 (x ,0)  =  0 


initial  condition 
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6(0, t) 
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boundary  conditions 


(25) 


where 


k(6) 


1  -  A0 


is  the  function  of  conductivity  with  temperature.  The  steps 
that  Luilcov  presents  to  find  the  analytic  solution  are 
lengthy  and  will  not  be  repeated  here.  The  solution  is 
given  by 


T  =  e/e, 


£  =  x/2 (aQT) 


a  =  X  a, 


3  =  -  In  (1  -  a) 


6  = 


C  = 


-u 


2  /  (u2-y  lr*  u2) "  2  du 

0 


(27) 


(28) 


(29) 


(30) 


(31) 


— 1— xt  C(u2-y  1  u2)'2-u]  expC/  (u2  -  y  lnu2)  ^  du,](32) 
(2„)^  n  0  1  1  1 


T  *  Higll  Cl-exp{-2  /  uj-y  lnuj)-15  du^D 

In  this  work,  ag  and  0g  are  always  unity  and  thus  0 
and  T  can  be  used  interchangeably. 


(33) 
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The  solution  of  these  equations  for  a  given  x  and 
t  is  numerical  and  is  not  trivial.  The  algorithm  used  by 
the  author  is  presented  in  the  Appendix. 
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III.  Method 


The  method  of  comparison  of  the  different  numerical 
methods  discussed  above  is  simple  in  principle.  All  that 
is  necessary  is  to  use  the  numerical  methods  on  Luikov's 
problem  and  compare  the  results  to  the  analytic  solution. 

The  peculiarities  of  the  test  problem  and  the  numerical 
methods,  however,  add  additional  considerations. 

The  problem  presented  by  Luikov  has  a  semi- inf inite 
domain.  This  is  always  difficult  to  treat  well  numerically, 
and  limitations  must  be  included  in  a  computer  code  that 
attempts  to  handle  a  non-finite  domain.  The  simplest 
approach  is  to  truncate  the  domain  at  some  point  and  halt 
the  calculations  when  that  area  of  the  domain  becomes  signi¬ 
ficant  to  the  problem.  This  usually  requires  subjective 
judgment  by  the  user.  For  Luikov's  problem,  where  the  driving 
force  is  the  heat  that  is  being  pumped  into  the  domain  at 
the  origin  (and  maintains  the  constant  temperature  at  that 
point),  it  is  reasonable  to  assume  that  temperatures  many 
orders  of  magnitude  smaller  than  those  at  the  origin  are 
probably  not  affecting  the  calculations  significantly.  This 
was  the  approach  used  by  the  author. 

Recall  that,  for  the  linear  heat  equation,  the  system 
matrices  for  finite  differences  and  finite  elements  are 
similar  in  many  cases  (Ref  9:416).  The  nonlinear  heat 
equation,  however,  yields  different  matrices  for  the  differ¬ 
ent  numerical  methods.  The  conductivity  in  Luikov's  problem 
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is  given  by 


a0 

kCe)  -  j-r0-^ 

where  0  >  A  >  1  .  The  author  set  =  1  for  all  test 

cases  because  it  is  just  a  normalization  factor.  As  A 
approaches  unity,  the  conductivity  becomes  very  sensitive 
to  temperature  and  therefore  highly  nonlinear.  To  emphasize 
differences  in  the  numerical  methods,  A  was  given  values 
relatively  close  to  unity. 

Another  difficulty  with  Luikov's  problem  is  the  dis¬ 
continuous  (in  time)  boundary  condition  at  the  origin. 

Others  have  found  this  to  cause  problems  with  numerical 
techniques  and  have  suggested  methods  to  overcome  it 
(Ref  1:18;  8:72).  This  author  believes  that  the  ability 
to  handle  this  type  of  boundary  condition  is  a  valid  test 
point  for  the  different  numerical  methods  and  has  included 
test  cases  that  use  it.  The  simplest  way  to  avoid  the 
problem  is  to  start  the  problem  at  a  later  time,  using  the 
analytic  solution  at  that  point  (Ref  8:72).  Test  cases  run 
in  this  manner,  also. 

The  major  parameters  that  can  be  varied  in  the  numerical 
solution  of  Luikov's  problem  are  the  time  step  length,  h  , 
the  size  to  the  spatial  increments,  Ax  ,  the  starting  time, 
and  A  .  The  latter  two  were  discussed  in  the  previous 
paragraphs.  The  different  methods  can  also  be  compared  at 
any  of  the  nodes  in  the  problem.  The  test  cases  were  created 


by  selecting  combinations  of  parameters  that  would  test 
various  qualities  of  the  numerical  methods.  The  author 
ran  many  other  cases  that  support  the  results  of  these  test 
cases,  but  the  data  from  those  are  not  presented  because  it 
would  be  voluminous  and  redundant.  The  parameters  for  each 
of  the  test  cases  are  shown  in  Table  1. 


TABLE 

1 

Test  Case  Parameters 

CASE 

PARAMETERS 

PURPOSE  OF  TEST 

1 

h  = 

4.0x10" 5 

Straight  solution  of 
Luikov's  problem  including 

Ax  = 

0.05 

the  discontinuous  boundary 

t0 

0.  0 

condition 

A  = 

0.90 

test 

node  =  2 

2 

h  = 

4.0x10" 5 

Same  as  first  test  except 

Ax  = 

0.05 

test  at  node  three 

t0  = 

0.0 

A  = 

0.  90 

test 

node  =  3 

3 

h  = 

2.0x10" 5 

Same  as  first  test  except 

Ax  = 

0.05 

A  =  0.95  and  thus  more 
nonlinear 

*0 

0.0 

A  = 

0.95 

test 

node  =  2 

4 

h  = 

4.0x10" 5 

Same  as  first  test  but 

Ax  = 

0.05 

started  at  t  >  0  with 

analytic  solution 

t0 

4.0x10 

A  = 

0.90 

test 

node  =  2 

5 

h  = 

2.0x10" 5 

Same  as  third  test  but 

Ax  = 

0.05 

started  at  t  >  0  with 

analytic  solution 

t0  = 

2.0x10 

A  « 

0.95 

test 

node  =  2 

time  step  length 
spatial  increment 
starting  time 


A  =  conductivity  parameter 
test  node  =  node  at  which  compar¬ 
ison  was  made 
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The  results  of  the  calculations  for  each  of  the  test 
cases  are  shown  in  Tables  2  through  6  and  Figures  2  through 
6.  The  first  test  case  includes  a  nonlinear  conductivity. 
The  conductivity  has  range  of  1.0  to  10.0,  depending  on  the 
normalized  temperature.  The  results  of  this  test  are 
significant.  The  lumped  mass  finite  element  method  followed 
the  analytic  solution  the  best  over  the  time  spectrum,  but 
the  distributed  mass  method  was  closer  at  t  =  4. Ox] O' 3 
The  finite  difference  methods  did  not  fare  as  well  as  the 
finite  element  methods. 


TEMPERATURE 


X  -  Analytic  0  -  Finite  Diff  3 

O-  Finite  Diff  1  Q  -  Finite  Elem  Distr  Mass 

A  -  Finite  Diff  2  +  -  Finite  Elem  Lumped  Mass 


Fig.  2.  Case  1  Results 


In  the  second  test  case,  the  third  node  was  chosen  as 
the  test  point.  It  was  thought  that  the  effect  of  the 
discontinuous  boundary  condition  would  be  less  at  this  point. 
The  lumped  mass  finite  element  method  was  the  best  performer 
over  the  whole  spectrum. 


s 

K 
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TABLE  3 


Finite  Element  1  =  Distributed  Mass 

2  =  Lumped  Mass 


TEMPERATURE 


Fig.  3.  Case  2  Results 
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In  the  third  case,  the  va]ue  of  A  was  increased  to 
0.95,  which  made  the  problem  more  nonlinear.  The  time  step 
length  was  decreased  accordingly  to  handle  the  higher  con- 
ducitvity.  All  the  methods  have  greater  difficulty  handling 
the  higher  conductivity,  but  the  limped  mass  finite  element 
method  still  is  the  best.  Finite  Difference  Scheme  1 
appears  to  do  well  at  early  times,  but  does  not  do  well  later 

-  3 

times  and  fails  to  converge  after  t  =  2.4  x  10 
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In  the  fourth  and  fifth  cases,  the  runs  were  repeated 
for  X  =  0.90  and  X  =  0.95  ,  except  the  time-stepping  was 

started  at  a  later  time  starting  with  the  exact  solution  at 
that  point.  The  behavior  in  these  cases  is  similar  to  those 
started  at  time  t  =  0  .  and  the  results  indicate  that  the 

lumped  mass  finite  element  method  does  the  best  over  the 
total  time  spectrum,  but  especially  at  later  times.  Note 
that  the  first  finite  difference  scheme  failed  to  converge 


TABLE 


Finite  Element  1  =  Distributed  Mass 

2  =  Lumped  Mass 
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Finite  Elem  Lumped  Mass 


Fig.  5.  Case  4  Results. 
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TEMPERATURE 


X-  Analytic  0-  Finite  Diff  3 

O-  Finite  Diff  1  ©-  Finite  Elem  Distr  Mass 

A-  Finite  Diff  2  +  -  Finite  Elem  Lumped  Mass 


Fig.  6.  Case  5  Results. 
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To  compare  the  speed  of  the  methods,  the  approximate 
time  it  took  the  methods  to  complete  ten  time  steps  was 
computed.  The  results  are  summarized  in  Table  7.  The  finite 
difference  schemes  were  the  fastest  and  the  lumped  mass 
finite  element  method  ran  a  close  third.  The  distributed 
mass  finite  difference  method  was  by  far  the  slowest.  It 
should  be  noted  that  an  anlaytic  form  of  the  integral  of 
the  conductivity  with  respect  to  temperature  was  used  in  the 
problem.  If  this  were  not  available,  it  is  possible  that 
the  finite  element  methods  would  be  even  slower. 

All  the  methods  used  a  Crank-Nicholson  scheme  for  the 
time  stepping.  Because  the  problem  is  nonlinear  and  Crank- 
Nicholson  is  implicit,  it  was  necessary  to  iterate  to  perform 
the  time  steps.  The  number  of  iterations  needed  to  perform 
a  time  step  was  kept  track  of  to  allow  comparisons  in  this 
area.  The  results  are  summarized  in  Table  7.  As  indicated 
by  this  table,  the  number  of  iterations  was  relatively  few 
for  the  methods,  with  the  distributed  mass  finite  element 
method  being  the  worst  case.  The  author  found  that  the  methods 
always  converged  except  for  first  finite  difference  scheme  in 
Cases  3  and  5. 
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TABLE  7 

Speed  and 

Convergence  Comparisons 

Time  to  Complete  Ten 

Iterations  per 

Method 

Time  Steps  (sec) 

Time 

Step 

Finite  Diff  1 

33 

2  - 

3 

Finite  Diff  2 

26 

2 

Finite  Diff  3 

30 

2 

Finite  Element 

Distributed  Mass 

43 

3  - 

4 

Finite  Element 

Lumped  Mass 

32 

2 

These  are  representative  numbers  based  on 
large  number  of  runs. 


V.  Conclusions 


The  results  from  the  test  cases  indicate  that  the  lumped 
mass  finite  element  method  is  far  and  away  the  best  choice 
from  among  those  tested  in  this  paper  in  terms  of  accuracy, 
although  the  distributed  mass  finite  element  method  was  some¬ 
times  better  at  late  times.  Because  the  speed  of  the  lumped 
mass  method  was  comparable  to  that  of  the  finite  difference 
methods,  there  is  no  basis  for  using  the  finite  difference 
methods  in  this  respect.  The  only  serious  problem  that  ould 
arise  is  if  the  integration  of  the  conductivity  with  respect 
to  temperature  is  difficult  or  time-consuming.  Because  it 
is  common  practice  to  use  Gaussian  quadrature  in  two-dimensional 
finite  elements,  however,  the  implementation  of  a  numerical 
integration  scheme  should  not  be  difficult.  From  these  results, 
it  is  concluded  that  the  lumped  mass  finite  element  method  is 
the  best  choice  from  among  those  tested  for  solving  the  non¬ 
linear  heat  equation. 
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|  VI.  Suggested  Extensions 

In  the  opinion  of  the  author,  the  conclusions  of  this 
work  are  striking.  The  fact  that  the  lumped  mass  finite 
l  element  method  proved  so  superior  was  not  anticipated  at  all 

This  casts  some  doubt  on  the  general  validity  of  the  work 
and  requires  other  trial  cases  to  be  examined.  Luikov  pre¬ 
sents  two  other  analytic  solutions  to  the  nonlinear  heat 
equation,  and  it  is  suggested  that  these  be  used  as  more 
benchmarks  for  the  methods  (Ref  7:502+,  504+) . 

Another  area  of  interest  would  be  the  comparison  of 


other  numerical  methods  to  those  tested  here.  Two  time¬ 
stepping  techniques  presented  by  Hughes  and  Trujillo  would 
be  of  particular  interest  to  this  problem  (Ref  6;  13). 


Bibliography 


1.  Blevins,  K.W.  A  Comparison  of  Finite  Difference  Methods 
for  the  SolutiorPoT  the  Transient  Heat  Conduct  ion  l-'quhtTon 
in  Inhomogeneous  Media.  Thesis.  Wright -Patterson  AFB  Oil: 
Air  Force  Institute  of  Technology,  March  1982. 

2.  Burden,  R.L.  ,  et_  al.  Numerical  Analysis.  Boston: 

Pr indie,  Weber  an3~Schmidt>  1978. 

3.  Clark,  M.R.  and  K.F.  Hansen.  Numerical  Methods  of 
Reactor  Analysis.  New  York:  Academic  Press,  1§64. 

4.  El-Wakil,  M.M.  Nuclear  Heat  Transfer.  LaGrange  Park, 

IL:  American  Nuclear  Society,  1971. 

5.  Fujita,  H.  "The  Exact  Pattern  of  a  Concentration- 
Dependent  Diffusion  in  Semi-Infinite  Medium,  Part  I." 
Textile  Research  Journal.  22 :  757-760  (November  1966). 

6.  Hughes,  T.J.R.  "Unconditionally  Stable  Algorithms  for 
Nonlinear  Heat  Conduction."  Computer  Methods  in  Applied 
Mechanics  and  Engineering.  10 :  135-139  (1977). 

7.  Luikov,  A.V.  Analytic  Heat  Diffusion  Theory.  New  York: 
Academic  Press,  1968. 

8.  Martin,  C.R.  An  Investigation  of  the  Numerical  Methods 
of  Finite  Differences  and  Finite  Elements  for  Digital 
Computer  Solution  of  the  Transient  Heat  Conduction 
(Diffusion)  Equation  Using  Optimum  Implicit  Formulations. 
Thesis .  Wright-Patterson  AFB  OH:  Air  Force  Institute  of 
Technology,  March  1978. 

9.  Myers,  G.E.  Analytical  Methods  in  Conduction  Heat 
Transfer .  New  York:  McGraw  Hill  Book  Co.,  1971. 

10.  Richtmeyer,  R.D.  and  K.W.  Morton.  Difference  Methods 
for  Initial-Value  Problems.  New  York :  Interscience 
Publishers,  1967. 

11.  Steam/Its  Generation  and  Use.  New  York:  Babcock  and 
Wilcox,  1978. 

12.  Trujillo,  D.M.  and  H.R.  Busby.  "Finite  Element  Nonlinear 

Heat  Transfer  Analysis  Using  a  Stable  Explicit  Method." 
Nuclear  Engineering  and  Design.  44 :  227-233  (1977). 

13.  Trujillo,  D.M.  "An  Unconditionally  Stable,  Explicit 

Algorithm  for  Finite  Element  Heat  Conduction  Analysis." 
Nuclear  Engineering  and  Design.  4 1 :  175-180  (1977). 

14.  Zienkiewicz,  O.C.  The  Finite  Element  Method.  London: 
McGraw  Hill  Book  Co.  (UK)  Limited,  1977. 


41 


APPENDIX  A 


Evaluation  of  the  Analytic  Solution 

To  compare  different  numerical  methods  on  their  ability 
to  solve  a  certain  class  of  problems,  it  is  often  necessary 
to  test  the  methods  on  a  sample  problem  from  the  class  that 
has  an  analytic  solution.  Unfortunately,  few  analytic 
solutions  are  known  for  nonlinear  heat  diffusion.  Luikov 
does  present  three  solutions  for  semi-infinite  domains. 

They  are  complicated  because  they  involve  numerical  integra¬ 
tion  and,  in  some  instances,  the  numerical  determination  of 
roots  of  functions  (Ref  7:499).  The  purpose  of  this  section 
is  to  demonstrate  how  to  implement  one  of  Luikov* s  solutions. 

The  time-dependent  nonlinear  heat  diffusion  equation  in 
one  dimension  is 


80 

3t 


& 


3jk 


(34) 


where  0  is  the  temperature,  x  is  the  time,  x  is  the 
spatial  coordinate,  and  k  is  the  temperature  dependent 
thermal  conductivity.  For  one  of  the  problems  presented 
by  Luikov,  the  conductivity  is  of  the  form 


k(9)  »  aQ  (1  -  X6) 


(35) 


where  X  and  ag  are  constants.  The  initial  and  boundary 
conditions  are 
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I  * 


e(x,o)  =  o 


6(0, t)  =  6. 


The  domain  of  this  problem  is  semi- inf inite  such  that  the 
left  hand  endpoint  is  at  the  origin  and  the  right  hand  end¬ 
point  is  at  infinity.  Through  numerous  variable  transforma¬ 
tions,  the  final  solution  is  given  as 


T  =  e/e, 


(37) 


t,  =  x/2  (aQx) 


(38) 


a  =  A  6, 


B  =  -ln(l-a) 


2  /  (u2-ylnu2) ~  2  du 
0 


(41) 


— — \-  [(u2-ylnu2)  ^-u]exp[/  [u^ - ylnu2 ]  ^Ddu.  (42) 

(2y)  2  0  1  1  1 

I',  exp  (■-  g)  Cl-exp{-2/  (u 2 -ylnu2 )  duj}]  (43) 


It  should  be  noted  that  all  these  equations  are  not  presented 
to  impress  the  reader,  but  actually  to  represent  the  most 
simplified  form  of  the  solution. 

Before  analyzing  these  equations,  it  is  important  to 
note  that  there  are  two  ways  to  approach  the  solution.  The 
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first  approach 

is 

to  find 

0  (or 

T  ) 

by  manipulating  the 

"solution"  above 

given  some 

spat  i; 

al  coordinate  x  and 

time  x  .  To 

do 

this,  one 

could 

use 

the  following  algorithm: 

Step 

1 

Find 

Z 

from 

Eq 

(38) 

Step 

2 

Find 

a 

from 

Eq 

(39) 

Step 

3 

Find 

3 

from 

Eq 

(40) 

Step 

4 

Find 

P 

from 

Eq 

(41) 

Step 

5 

Find 

u 

from 

Eq 

(42) 

Step 

6 

Find 

T 

from 

Eq 

(43) 

Step 

7 

Find 

e 

from 

Eq 

(37) 

The  first  three  steps  of  this  algorithm  are  rather  trivial. 

Steps  4,  5,  and  6,  however,  require  numerical  integration 
and  Steps  4  and  5  also  require  numerical  root  finding. 

Obviously,  a  price  is  paid  for  wanting  the  solution  at  a 
particular  point.  The  second  approach  is  to  select  various 
values  of  u  and  y  and  backtrack  to  find  the  corresponding 
values  of  x  ,  t  ,  3  and  a  .  With  judicious  selection 
of  u  and  y  ,  one  can  get  a  good  idea  of  the  behavior  of 
the  solution,  although  it  would  take  a  stroke  of  good  luck 
to  find  values  at  predetermined  values  of  x  ,  x  and  a 
This  second  approach  would  still  require  numerical  integration, 
but  would  leave  out  the  numerical  root  finding.  For  the  pur¬ 
poses  of  computer  code  comparison,  however,  the  first  approach 
is  more  appropriate.  It  will  be  demonstrated  that  the  numerical 
integration  requires  much  more  sophistication  and  computer 
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time  than  the  numerical  root  finding  anyway,  and  that  the 
root  finding  is  fairly  straightforward. 

The  integral  that  must  be  evaluated  is  the  equation 


/  (u2  -  y)?nu2)  ^  du 

0 


Note  the  integrands  in  Eqs  (41),  (42),  and  (43)  are  identical, 
which  will  be  shown  to  be  greatly  advantageous  below.  Note 
that  there  is  a  logarithm  in  the  integrand.  Although  the 
integrand  itself  behaves  nicely  as  it  approaches  zero,  the 
calculation  of  the  logarithm  becomes  difficult.  The  first 
task,  then,  is  to  determine  how  close  to  the  origin  the  left 
hand  limit  of  the  integration  must  be  to  limit  the  error  to 
some  reasonable  level.  The  integrals  in  Eqs  (41),  (42),  and 
(43)  can  be  rewritten 


=  /  [u2-yln(u2)3  2  dux  +  /  [u^-yln  (u^ 2] 


2~\~h 


Define  I  as 
e 


,  min  -v 

Ie  “  /  [u2-yln(uj  )  3  2  dux 


Thus,  I  is  an  absolute  measure  of  the  error  that  results 
e 

from  fixing  the  left  hand  limit  at  um^n  instead  of  the 
origin.  Now 
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hi  • 

f  ' 

* 

— ■ — ■ — — ■"* 

y  s  0  ,  0  <  u  <  1 

(47) 

1  * 

implies 

i . 

f 

so 

-  y  In  (u  2  )  >  0 

(48) 

> 

Cu2 

-ylnu2U^  >  C-ylnu2!!"2 

(49) 

*  < 

Cu2 

-ylnu23  ^  <  [-ylnu2D 

(50) 

hence 

I 

e 

u  .  u  . 

mm  ,  mm  i 

=  /  [u2 -ylnu2]"^  du,  <  /  C-ylnu2]"2  du 

0  11  1  0 

(51) 

1  1 

»_ ' 

I 

e 

,  Umin  , 

<  —  /  C-lnur*  du 

ST\x  0 

(52) 

►  — — 

If 

0  <,  u  <  u  . 

mm 

(53) 

<  w 

*  • 

then 

i  ,* 

\ 

• 

-In 

u  <  In  u  • 
mm 

(54) 

!»* 

•*,  • 

1  ♦ 
r  v 
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and 


(55) 


/-inu  /^lnu 

mm 


Thus 


u  . 
mm 


J1I-LJ1  -t 

I  7---  du  <  f 

0  /- lnu  0 


u  . 
mm 


/^inu  . 
mm 


(56) 


-lnu 


u  . 
mm 


I  <  — 

e  /2TT 


umin 

"lnumin 


u  ,  <  1  ,  /-lnu  .  >  1 

mm  ’  mm 


(59) 


one  gets 


u  . 

X  <  min 

e  /2TT 


(60) 


Suppose 


I(u)  >  lx  10' 
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For  five  significant  figures  of  accuracy, 


u 


mm 

/2p 


1x10' 5  I(u) 


<  1x10 


-  8 


(62  ) 


and 


u  . 

mm  < 

,/2V 


1x10 


(63) 


If  y  is  bounded  by 


y  >  0.001 


(64) 


then 


4.47  x  10~10 


u. 


min 


(65) 


or 


u  .  <  1  x  10" 1 °  (66) 

mm 

Hence  a  value  of  u  .  has  been  determined  that  will  limit 

mm 

the  error  to  five  significant  figures.  (Actually,  most  of 
the  bounds  on  different  variables  were  very  conservative  and 
the  error  should  be  much  less  than  this.) 

Now  that  practical  limits  of  integration  have  been 
found,  it  is  necessary  to  select  a  numerical  integrator  to 
perform  the  integration.  The  Romberg  method  as  described  by 


1  , 


B 


t 


r. 

P 


i 


Burden  was  used  by  the  author  (Ref  2:206+).  The  method 
provides  an  indication  of  the  error  of  the  approximation  it 
finds  for  the  integral.  To  get  an  accuracy  of  about  five 
significant  figures,  the  author  found  it  necessary  to  make 
as  many  as  16,385  function  evaluations  over  the  limits  from 
x  =  10”10  to  x  =  1.0  .  This  many  evaluations  requires  a 
considerable  amount  of  computer  time  in  itself,  but  in  addition, 
the  program  had  to  be  executed  in  double  precision  (i.e., 
microcomputer  double  precision,  which  is  about  fifteen  signi¬ 
ficant  figures)  which  slowed  execution  further.  The  inte¬ 
grator  appeared  to  have  the  most  trouble  near  the  origin. 

More  will  be  said  about  this  later. 

Recall  the  goal  at  this  point  is  to  implement  Step  4 
of  the  algorithm  above,  i.e.,  find  y  using  Eq  (41)  given 
some  3  .  This  can  be  done  with  little  trouble  using  the 

Newton-Raphson  Method  (Ref  2:39+).  Equation  (41)  can  be 
rewritten 


F(y)  =  2 J  Cu2-plnu23  du  -  8 


(67) 


where  one  wants  to  find  y  such  that  F(y)  =  0  .  Using 
the  standard  Newton-Raphson  Method,  an  iterative  equation 
to  determine  y  is 


'i»l  *  “i  '  FCm^/D^O..) 


(68) 


f 
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The  derivative  is  just 


D  F  (y)  =  /  (u2-  lnu2)  3/2(lnu2)  du  (69) 

p  0 

It  turns  out  that  this  integral  converges  more  quickly  than 
the  previous  one.  The  Newton-Raphson  Method  also  converges 
nicely  to  better  than  five  significant  figures  within  five 
or  six  iterations.  The  problem  is  that  two  integrals  have 
to  be  evaluated  for  each  iteration.  Note  that  X  ,  and 
hence  a  and  3  ,  are  really  parameters  for  a  particular 

case  and  not  variables.  For  a  particular  choice  of  X  , 
one  has  to  evaluate  y  only  once. 

Unfortunately,  one  would  have  to  execute  Steps  5  and  6 
of  the  algorithm  many  times  to  get  an  idea  of  what  the  solu¬ 
tion  to  the  problem  really  looks  like  and  to  compare  numerical 
methods  at  more  than  one  point  in  time  and  space.  Step  5  is 
similar  to  Step  4  because  of  the  combination  of  numerical 
root  finding  and  integral  evaluation  that  must  be  used.  The 
integral  evaluations,  however,  can  be  made  much  easier  in 
these  steps  than  in  the  previous  one.  Recall  that  the  trouble 
spot  in  the  numerical  integration  is  near  the  origin.  The 
integral  to  be  evaluated  can  be  rewritten 

/  (u? -ylnu? ) "^du,  =  /  (u2 -ylmi 2 ) '**  du.. 

0  x  0  1  1  1 

1  . 

+  -/  (u2-ylnu2  )"*  du,  (70) 

u  1  1  1 
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Using  Eq  (41),  one  can  substitute  for  the  first  term  on 
the  right  hand  side  to  get 

/  (u^-ylnup”'1  duj  =  -j~  -  f  (u*-ylnup’'s  du1  (71) 

Note  that  no  approximations  have  been  made  because  the  exact 
value  of  B  can  be  found  for  all  intents  and  purposes.  It 
turns  out  that  this  integral  is  much  easier  to  evaluate  for 
the  same  accuracy  than  the  previous  one  because  it  stays  away 
from  the  origin.  Finally,  the  Newton-Raphson  Method  performs 
just  as  nicely  on  Eq  (42)  as  it  did  on  Eq  (41). 

With  the  simplifications  that  have  been  suggested,  the 
final  algorithm  is  expanded  as  shown  on  the  following  page. 

Table  8  presents  various  values  of  a  ,  3  , 


and  y 


Step 

1 

Find 

a 

from  Eq  (39) 

Step 

2 

Find 

8 

from  Eq  (40) 

Step 

3 

Find 

V 

from  Newton-Raphson  and  Romberg  Methods  on 

B 

=  2 

1 

/ 

(u2  -  plnu2)  1  du 

1x10*  1  ° 


Step  4  For  every  pair  of  (x,x),  do 
4a  Find  £  from  Eq  (38) 

4b  Find  u  by  using  Newton-Raphson  and  Romberg  on 
5  -  — C(u2  -  vlrm2)h  -  u3 

(2|i  )* 

x  exp[|-  -  /  (u2  -  plnu2)"2  duj] 

4c  Find  T  by  using  Romberg  on 

T  =  1 

1  -  exp  (-6) 

0  1  1 

x  (1  -  exp[-2(^  -  /  (u2  -  plnu2)  *  du^)]} 

4d  Find  0  from  Eq  (37) 


Algorithm  to  Numerically  Evaluate 
the  Analytic  Solution  for  an  Example 
of  the  Nonlinear  Heat  Equation 
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TABLE  8 


Values  of  a  , 

8  ,  and  y 

for  the  Test  Problem 

a 

8 

y 

0.99 

4.605165 

0.009684546 

0.98 

3.911202 

0.02360335 

0.975 

3.688878 

0.03174528 

0.95 

2.995732 

0.08306863 

0.90 

2.302585 

0.2385708 

0.80 

1.609438 

0.8125746 
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